### heterogeneous treatment effects left/right ###


## create three-way lr dummi with moderates

# create lr dummi and subdata
ESS10$lmrdummi[ESS10$lrscale <5] <- "left"
ESS10$lmrdummi[ESS10$lrscale >5] <- "right"
ESS10$lmrdummi[ESS10$lrscale == 5] <- "moderate"
table(ESS10$lmrdummi)

left10 <- ESS10 %>% 
  filter(ESS10$lmrdummi == "left")
right10 <- ESS10 %>% 
  filter(ESS10$lmrdummi == "right")
moderate10 <- ESS10 %>% 
  filter(ESS10$lmrdummi == "moderate")

migleft <- RDestimate(mig ~ date | cntry, data = left10, cutpoint = 0, bw = 14)
summary(migleft)
lowmigleft <- migleft$ci[1,1]
upmigleft <- migleft$ci[1,2]
coefmigleft <- migleft$est[1]

migmod <- RDestimate(mig ~ date | cntry, data = moderate10, cutpoint = 0, bw = 14)
summary(migmod)
lowmigmod <- migmod$ci[1,1]
upmigmod <- migmod$ci[1,2]
coefmigmod <- migmod$est[1]

migright <- RDestimate(mig ~ date | cntry, data = right10, cutpoint = 0, bw = 14)
summary(migright)
lowmigright <- migright$ci[1,1]
upmigright <- migright$ci[1,2]
coefmigright <- migright$est[1]

ecoleft <- RDestimate(eco ~ date | cntry, data = left10, cutpoint = 0, bw = 14)
summary(ecoleft)
lowecoleft <- ecoleft$ci[1,1]
upecoleft <- ecoleft$ci[1,2]
coefecoleft <- ecoleft$est[1]

ecomod <- RDestimate(eco ~ date | cntry, data = moderate10, cutpoint = 0, bw = 14)
summary(ecomod)
lowecomod <- ecomod$ci[1,1]
upecomod <- ecomod$ci[1,2]
coefecomod <- ecomod$est[1]

ecoright <- RDestimate(eco ~ date | cntry, data = right10, cutpoint = 0, bw = 14)
summary(ecoright)
lowecoright <- ecoright$ci[1,1]
upecoright <- ecoright$ci[1,2]
coefecoright <- ecoright$est[1]

culleft <- RDestimate(cul ~ date | cntry, data = left10, cutpoint = 0, bw = 14)
summary(culleft)
lowculleft <- culleft$ci[1,1]
upculleft <- culleft$ci[1,2]
coefculleft <- culleft$est[1]

culmod <- RDestimate(cul ~ date | cntry, data = moderate10, cutpoint = 0, bw = 14)
summary(culmod)
lowculmod <- culmod$ci[1,1]
upculmod <- culmod$ci[1,2]
coefculmod <- culmod$est[1]

culright <- RDestimate(cul ~ date | cntry, data = right10, cutpoint = 0, bw = 14)
summary(culright)
lowculright <- culright$ci[1,1]
upculright <-culright$ci[1,2]
coefculright <- culright$est[1]

hetcoef <- c(coefmigleft, coefmigmod, coefmigright, coefecoleft, coefecomod, 
             coefecoright, coefculleft, coefculmod, coefculright)
hetcoef <- as.data.frame(hetcoef)

hetup <- c(upmigleft, upmigmod, upmigright, upecoleft, upecomod, upecoright, 
           upculleft, upculmod, upculright)
hetup <- as.data.frame(hetup)

hetlow <- c(lowmigleft, lowmigmod, lowmigright, lowecoleft, lowecomod, lowecoright, 
            lowculleft, lowculmod, lowculright)
hetlow <- as.data.frame(hetlow)

hetplot <- bind_cols(hetcoef,hetlow,hetup)
hetplot$Politically <- as.character(c("Leftwing", "Moderate", "Rightwing", "Leftwing", "Moderate", "Rightwing", "Leftwing", "Moderate", "Rightwing"))
hetplot$var <- as.character(c("Overall", "Overall", "Overall", "Economy", "Economy", "Economy", "Culture", "Culture", "Culture"))
level_order <- c("Overall", "Economy", "Culture")

leftmodright <- ggplot(hetplot, aes(x = var, y = hetcoef, color = Politically)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 6, position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin = hetlow, ymax = hetup), 
                size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-3,3) +
  theme_light(base_size = 30) +
  xlab("") +
  ylab("") +
  scale_color_manual(values = c("#F8766D","#00BA38", "#00A9FF")) +
  scale_x_discrete(limits = level_order)

leftmodright




### heterogeneous country effects ### 
table(smol$cntry)

bg <- ESS10 %>% # had troops
  filter(cntry == "BG")

ch <- ESS10 %>% # had troops 
  filter(cntry == "CH")

cz <- ESS10 %>% # had troops 
  filter(cntry == "CZ")

ee <- ESS10 %>% # had troops 
  filter(cntry == "EE")

fr <- ESS10 %>% # had troops
  filter(cntry == "FR")

hr <- ESS10 %>% # had troops
  filter(cntry == "HR")

hu <- ESS10 %>% # had troops
  filter(cntry == "HU")

is <- ESS10 %>% # had troops 
  filter(cntry == "IS")

lt <- ESS10 %>% # had troops
  filter(cntry == "LT")

no <- ESS10 %>% # had troops 
  filter(cntry == "NO")

pt <- ESS10 %>% # had troops 
  filter(cntry == "PT")

si <- ESS10 %>% # had troops
  filter(cntry == "SI")

sk <- ESS10 %>% # had troops
  filter(cntry == "SK")


## bulgaria
rd_migbg <- RDestimate(mig ~ date, data = bg, cutpoint = 0, bw = 14)
summary(rd_migbg)
rd_ecobg <- RDestimate(eco ~ date, data = bg, cutpoint = 0, bw = 14)
summary(rd_ecobg)
rd_culbg <- RDestimate(cul ~ date, data = bg, cutpoint = 0, bw = 14)
summary(rd_culbg)

lowmigbg <- rd_migbg$ci[1,1]
upmigbg <- rd_migbg$ci[1,2]
coefmigbg <- rd_migbg$est[1]

lowecobg <- rd_ecobg$ci[1,1]
upecobg <- rd_ecobg$ci[1,2]
coefecobg <- rd_ecobg$est[1]

lowculbg <- rd_culbg$ci[1,1]
upculbg <- rd_culbg$ci[1,2]
coefculbg <- rd_culbg$est[1]


## switzerland
rd_migch <- RDestimate(mig ~ date, data = ch, cutpoint = 0, bw = 14)
summary(rd_migch)
rd_ecoch <- RDestimate(eco ~ date, data = ch, cutpoint = 0, bw = 14)
summary(rd_ecoch)
rd_culch <- RDestimate(cul ~ date, data = ch, cutpoint = 0, bw = 14)
summary(rd_culch)

lowmigch <- rd_migch$ci[1,1]
upmigch <- rd_migch$ci[1,2]
coefmigch <- rd_migch$est[1]

lowecoch <- rd_ecoch$ci[1,1]
upecoch <- rd_ecoch$ci[1,2]
coefecoch <- rd_ecoch$est[1]

lowculch <- rd_culch$ci[1,1]
upculch <- rd_culch$ci[1,2]
coefculch <- rd_culch$est[1]


## czech republic
rd_migcz <- RDestimate(mig ~ date, data = cz, cutpoint = 0, bw = 14)
summary(rd_migcz)
rd_ecocz <- RDestimate(eco ~ date, data = cz, cutpoint = 0, bw = 14)
summary(rd_ecocz)
rd_culcz <- RDestimate(cul ~ date, data = cz, cutpoint = 0, bw = 14)
summary(rd_culcz)

lowmigcz <- rd_migcz$ci[1,1]
upmigcz <- rd_migcz$ci[1,2]
coefmigcz <- rd_migcz$est[1]

lowecocz <- rd_ecocz$ci[1,1]
upecocz <- rd_ecocz$ci[1,2]
coefecocz <- rd_ecocz$est[1]

lowculcz <- rd_culcz$ci[1,1]
upculcz <- rd_culcz$ci[1,2]
coefculcz <- rd_culcz$est[1]


## estonia
rd_migee <- RDestimate(mig ~ date, data = ee, cutpoint = 0, bw = 14)
summary(rd_migee)
rd_ecoee <- RDestimate(eco ~ date, data = ee, cutpoint = 0, bw = 14)
summary(rd_ecoee)
rd_culee <- RDestimate(cul ~ date, data = ee, cutpoint = 0, bw = 14)
summary(rd_culee)

lowmigee <- rd_migee$ci[1,1]
upmigee <- rd_migee$ci[1,2]
coefmigee <- rd_migee$est[1]

lowecoee <- rd_ecoee$ci[1,1]
upecoee <- rd_ecoee$ci[1,2]
coefecoee <- rd_ecoee$est[1]

lowculee <- rd_culee$ci[1,1]
upculee <- rd_culee$ci[1,2]
coefculee <- rd_culee$est[1]


## france
# sample too small for RDD (34 cases at 14 day bw)


## croatia
rd_mighr <- RDestimate(mig ~ date, data = hr, cutpoint = 0, bw = 14)
summary(rd_mighr)
rd_ecohr <- RDestimate(eco ~ date, data = hr, cutpoint = 0, bw = 14)
summary(rd_ecohr)
rd_culhr <- RDestimate(cul ~ date, data = hr, cutpoint = 0, bw = 14)
summary(rd_culhr)

lowmighr <- rd_mighr$ci[1,1]
upmighr <- rd_mighr$ci[1,2]
coefmighr <- rd_mighr$est[1]

lowecohr <- rd_ecohr$ci[1,1]
upecohr <- rd_ecohr$ci[1,2]
coefecohr <- rd_ecohr$est[1]

lowculhr <- rd_culhr$ci[1,1]
upculhr <- rd_culhr$ci[1,2]
coefculhr <- rd_culhr$est[1]


## hungary
rd_mighu <- RDestimate(mig ~ date, data = hu, cutpoint = 0, bw = 14)
summary(rd_mighu)
rd_ecohu <- RDestimate(eco ~ date, data = hu, cutpoint = 0, bw = 14)
summary(rd_ecohu)
rd_culhu <- RDestimate(cul ~ date, data = hu, cutpoint = 0, bw = 14)
summary(rd_culhu)

lowmighu <- rd_mighu$ci[1,1]
upmighu <- rd_mighu$ci[1,2]
coefmighu <- rd_mighu$est[1]

lowecohu <- rd_ecohu$ci[1,1]
upecohu <- rd_ecohu$ci[1,2]
coefecohu <- rd_ecohu$est[1]

lowculhu <- rd_culhu$ci[1,1]
upculhu <- rd_culhu$ci[1,2]
coefculhu <- rd_culhu$est[1]


## iceland
rd_migis <- RDestimate(mig ~ date, data = is, cutpoint = 0, bw = 14)
summary(rd_migis)
rd_ecois <- RDestimate(eco ~ date, data = is, cutpoint = 0, bw = 14)
summary(rd_ecois)
rd_culis <- RDestimate(cul ~ date, data = is, cutpoint = 0, bw = 14)
summary(rd_culis)

lowmigis <- rd_migis$ci[1,1]
upmigis <- rd_migis$ci[1,2]
coefmigis <- rd_migis$est[1]

lowecois <- rd_ecois$ci[1,1]
upecois <- rd_ecois$ci[1,2]
coefecois <- rd_ecois$est[1]

lowculis <- rd_culis$ci[1,1]
upculis <- rd_culis$ci[1,2]
coefculis <- rd_culis$est[1]


## lithuania
rd_miglt <- RDestimate(mig ~ date, data = lt, cutpoint = 0, bw = 14)
summary(rd_miglt)
rd_ecolt <- RDestimate(eco ~ date, data = lt, cutpoint = 0, bw = 14)
summary(rd_ecolt)
rd_cullt <- RDestimate(cul ~ date, data = lt, cutpoint = 0, bw = 14)
summary(rd_cullt)

lowmiglt <- rd_miglt$ci[1,1]
upmiglt <- rd_miglt$ci[1,2]
coefmiglt <- rd_miglt$est[1]

lowecolt <- rd_ecolt$ci[1,1]
upecolt <- rd_ecolt$ci[1,2]
coefecolt <- rd_ecolt$est[1]

lowcullt <- rd_cullt$ci[1,1]
upcullt <- rd_cullt$ci[1,2]
coefcullt <- rd_cullt$est[1]

## norway 
rd_migno <- RDestimate(mig ~ date, data = no, cutpoint = 0, bw = 14)
summary(rd_migno)
rd_econo <- RDestimate(eco ~ date, data = no, cutpoint = 0, bw = 14)
summary(rd_econo)
rd_culno <- RDestimate(cul ~ date, data = no, cutpoint = 0, bw = 14)
summary(rd_culno)

lowmigno <- rd_migno$ci[1,1]
upmigno <- rd_migno$ci[1,2]
coefmigno <- rd_migno$est[1]

lowecono <- rd_econo$ci[1,1]
upecono <- rd_econo$ci[1,2]
coefecono <- rd_econo$est[1]

lowculno <- rd_culno$ci[1,1]
upculno <- rd_culno$ci[1,2]
coefculno <- rd_culno$est[1]

## portugal
# sample too small for RDD (59 cases at 14 day bw)


## slovenia
# sample too small for RDD (21 cases at 14 day bw)


## slovakia
rd_migsk <- RDestimate(mig ~ date, data = sk, cutpoint = 0, bw = 14)
summary(rd_migsk)
rd_ecosk <- RDestimate(eco ~ date, data = sk, cutpoint = 0, bw = 14)
summary(rd_ecosk)
rd_culsk <- RDestimate(cul ~ date, data = sk, cutpoint = 0, bw = 14)
summary(rd_culsk)

lowmigsk <- rd_migsk$ci[1,1]
upmigsk <- rd_migsk$ci[1,2]
coefmigsk <- rd_migsk$est[1]

lowecosk <- rd_ecosk$ci[1,1]
upecosk <- rd_ecosk$ci[1,2]
coefecosk <- rd_ecosk$est[1]

lowculsk <- rd_culsk$ci[1,1]
upculsk <- rd_culsk$ci[1,2]
coefculsk <- rd_culsk$est[1]



ccoef <- c(coefmigbg, coefecobg, coefculbg,
           coefmigch, coefecoch, coefculch,
           coefmigcz, coefecocz, coefculcz,
           coefmigee, coefecoee, coefculee,
           coefmighr, coefecohr, coefculhr,
           coefmighu, coefecohu, coefculhu,
           coefmigis, coefecois, coefculis,
           coefmiglt, coefecolt, coefcullt,
           coefmigno, coefecono, coefculno,
           coefmigsk, coefecosk, coefculsk)
ccoef <- as.data.frame(ccoef)

cup <- c(upmigbg, upecobg, upculbg,
         upmigch, upecoch, upculch,
         upmigcz, upecocz, upculcz,
         upmigee, upecoee, upculee,
         upmighr, upecohr, upculhr,
         upmighu, upecohu, upculhu,
         upmigis, upecois, upculis,
         upmiglt, upecolt, upcullt,
         upmigno, upecono, upculno,
         upmigsk, upecosk, upculsk)
cup <- as.data.frame(cup)

clow <- c(lowmigbg, lowecobg, lowculbg,
          lowmigch, lowecoch, lowculch,
          lowmigcz, lowecocz, lowculcz,
          lowmigee, lowecoee, lowculee,
          lowmighr, lowecohr, lowculhr,
          lowmighu, lowecohu, lowculhu,
          lowmigis, lowecois, lowculis,
          lowmiglt, lowecolt, lowcullt,
          lowmigno, lowecono, lowculno,
          lowmigsk, lowecosk, lowculsk)
clow <- as.data.frame(clow)

cplot <- bind_cols(ccoef,clow,cup)
cplot$Country <- as.character(c("Bulgaria", "Bulgaria", "Bulgaria",
                                "Switzerland", "Switzerland", "Switzerland",
                                "Czech Republic", "Czech Republic", "Czech Republic",
                                "Estonia", "Estonia", "Estonia",
                                "Croatia", "Croatia", "Croatia",
                                "Hungary", "Hungary", "Hungary",
                                "Iceland", "Iceland", "Iceland",
                                "Lithuania", "Lithuania", "Lithuania",
                                "Norway", "Norway", "Norway",
                                "Slovakia", "Slovakia", "Slovakia"))

cplot$Outcome <- as.character(c("Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture",
                                "Overall", "Economy", "Culture"))

cplot$Country <- factor(cplot$Country)

cplot$Outcome <- factor(cplot$Outcome)
cplot$Outcome <- factor(cplot$Outcome, levels = rev(levels(cplot$Outcome)))


cntryplot <- ggplot(cplot, aes(x = Country, y = ccoef, color = Outcome)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 6, position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin = clow, ymax = cup), 
                size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-5,5) +
  theme_light(base_size = 25) +
  xlab("Country") +
  ylab("RDD Coefficient") +
  scale_color_manual(values = c("#00BA38", "#00A9FF", "#F8766D")) 
cntryplot


# heterogeneous region effects

ESS10$eastwest <- NA
ESS10$eastwest[ESS10$cntry == "BG" | ESS10$cntry == "CZ" |
                 ESS10$cntry == "EE" | ESS10$cntry == "HR" |
                 ESS10$cntry == "HU" | ESS10$cntry == "LT" |
                 ESS10$cntry == "SI" | ESS10$cntry == "SK"] <- "East"

ESS10$eastwest[ESS10$cntry == "CH" | ESS10$cntry == "FR" |
                 ESS10$cntry == "IS" | ESS10$cntry == "NO" |
                 ESS10$cntry == "PT"] <- "West"


east10 <- ESS10 %>% 
  filter(ESS10$eastwest == "East")
west10 <- ESS10 %>% 
  filter(ESS10$eastwest == "West")

migeast <- RDestimate(mig ~ date | cntry, data = east10, cutpoint = 0, bw = 14)
summary(migeast)
lowmigeast <- migeast$ci[1,1]
upmigeast <- migeast$ci[1,2]
coefmigeast <- migeast$est[1]

migwest <- RDestimate(mig ~ date | cntry, data = west10, cutpoint = 0, bw = 14)
summary(migwest)
lowmigwest <- migwest$ci[1,1]
upmigwest <- migwest$ci[1,2]
coefmigwest <- migwest$est[1]

ecoeast <- RDestimate(eco ~ date | cntry, data = east10, cutpoint = 0, bw = 14)
summary(ecoeast)
lowecoeast <- ecoeast$ci[1,1]
upecoeast <- ecoeast$ci[1,2]
coefecoeast <- ecoeast$est[1]

ecowest <- RDestimate(eco ~ date | cntry, data = west10, cutpoint = 0, bw = 14)
summary(ecowest)
lowecowest <- ecowest$ci[1,1]
upecowest <- ecowest$ci[1,2]
coefecowest <- ecowest$est[1]

culeast <- RDestimate(cul ~ date | cntry, data = east10, cutpoint = 0, bw = 14)
summary(culeast)
lowculeast <- culeast$ci[1,1]
upculeast <- culeast$ci[1,2]
coefculeast <- culeast$est[1]

culwest <- RDestimate(cul ~ date | cntry, data = west10, cutpoint = 0, bw = 14)
summary(culwest)
lowculwest <- culwest$ci[1,1]
upculwest <-culwest$ci[1,2]
coefculwest <- culwest$est[1]

hetcoefr <- c(coefmigeast, coefmigwest, coefecoeast, 
              coefecowest, coefculeast, coefculwest)
hetcoefr <- as.data.frame(hetcoefr)

hetupr <- c(upmigeast, upmigwest, upecoeast, upecowest, 
            upculeast, upculwest)
hetupr <- as.data.frame(hetupr)

hetlowr <- c(lowmigeast, lowmigwest, lowecoeast, lowecowest, 
             lowculeast, lowculwest)
hetlowr <- as.data.frame(hetlowr)

hetplotr <- bind_cols(hetcoefr,hetlowr,hetupr)
hetplotr$Region <- as.character(c("East", "West", "East", "West", "East", "West"))
hetplotr$varr <- as.character(c("Overall", "Overall", "Economy", "Economy", "Culture", "Culture"))
level_order <- c("Overall", "Economy", "Culture")

region <- ggplot(hetplotr, aes(x = varr, y = hetcoefr, color = Region)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 6, position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin = hetlowr, ymax = hetupr), 
                size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-3,3) +
  xlab("") +
  ylab("") +
  theme_light(base_size = 30) +
  scale_color_manual(values = c("#FB61D7","#00BE67")) +
  scale_x_discrete(limits = level_order)
region

ggarrange(leftmodright, region)



# by region and ideology combined

ESS10$both <- NA
ESS10$both[ESS10$eastwest == "East" & ESS10$lmrdummi == "left"] <- "Eastern Left"
ESS10$both[ESS10$eastwest == "East" & ESS10$lmrdummi == "moderate"] <- "Eastern Moderate"
ESS10$both[ESS10$eastwest == "East" & ESS10$lmrdummi == "right"] <- "Eastern Right"
ESS10$both[ESS10$eastwest == "West" & ESS10$lmrdummi == "left"] <- "Western Left"
ESS10$both[ESS10$eastwest == "West" & ESS10$lmrdummi == "moderate"] <- "Western Moderate"
ESS10$both[ESS10$eastwest == "West" & ESS10$lmrdummi == "right"] <- "Western Right"

eastl10 <- ESS10 %>% 
  filter(ESS10$both == "Eastern Left")
eastm10 <- ESS10 %>% 
  filter(ESS10$both == "Eastern Moderate")
eastr10 <- ESS10 %>% 
  filter(ESS10$both == "Eastern Right")

westl10 <- ESS10 %>% 
  filter(ESS10$both == "Western Left")
westm10 <- ESS10 %>% 
  filter(ESS10$both == "Western Moderate")
westr10 <- ESS10 %>% 
  filter(ESS10$both == "Western Right")


migeastl <- RDestimate(mig ~ date | cntry, data = eastl10, cutpoint = 0, bw = 14)
summary(migeastl)
lowmigeastl <- migeastl$ci[1,1]
upmigeastl <- migeastl$ci[1,2]
coefmigeastl <- migeastl$est[1]

migwestl <- RDestimate(mig ~ date | cntry, data = westl10, cutpoint = 0, bw = 14)
summary(migwestl)
lowmigwestl <- migwestl$ci[1,1]
upmigwestl <- migwestl$ci[1,2]
coefmigwestl <- migwestl$est[1]

ecoeastl <- RDestimate(eco ~ date | cntry, data = eastl10, cutpoint = 0, bw = 14)
summary(ecoeastl)
lowecoeastl <- ecoeastl$ci[1,1]
upecoeastl <- ecoeastl$ci[1,2]
coefecoeastl <- ecoeastl$est[1]

ecowestl <- RDestimate(eco ~ date | cntry, data = westl10, cutpoint = 0, bw = 14)
summary(ecowestl)
lowecowestl <- ecowestl$ci[1,1]
upecowestl <- ecowestl$ci[1,2]
coefecowestl <- ecowestl$est[1]

culeastl <- RDestimate(cul ~ date | cntry, data = eastl10, cutpoint = 0, bw = 14)
summary(culeastl)
lowculeastl <- culeastl$ci[1,1]
upculeastl <- culeastl$ci[1,2]
coefculeastl <- culeastl$est[1]

culwestl <- RDestimate(cul ~ date | cntry, data = westl10, cutpoint = 0, bw = 14)
summary(culwestl)
lowculwestl <- culwestl$ci[1,1]
upculwestl <-culwestl$ci[1,2]
coefculwestl <- culwestl$est[1]



migeastm <- RDestimate(mig ~ date | cntry, data = eastm10, cutpoint = 0, bw = 14)
summary(migeastm)
lowmigeastm <- migeastm$ci[1,1]
upmigeastm <- migeastm$ci[1,2]
coefmigeastm <- migeastm$est[1]

migwestm <- RDestimate(mig ~ date | cntry, data = westm10, cutpoint = 0, bw = 14)
summary(migwestm)
lowmigwestm <- migwestm$ci[1,1]
upmigwestm <- migwestm$ci[1,2]
coefmigwestm <- migwestm$est[1]

ecoeastm <- RDestimate(eco ~ date | cntry, data = eastm10, cutpoint = 0, bw = 14)
summary(ecoeastm)
lowecoeastm <- ecoeastm$ci[1,1]
upecoeastm <- ecoeastm$ci[1,2]
coefecoeastm <- ecoeastm$est[1]

ecowestm <- RDestimate(eco ~ date | cntry, data = westm10, cutpoint = 0, bw = 14)
summary(ecowestm)
lowecowestm <- ecowestm$ci[1,1]
upecowestm <- ecowestm$ci[1,2]
coefecowestm <- ecowestm$est[1]

culeastm <- RDestimate(cul ~ date | cntry, data = eastm10, cutpoint = 0, bw = 14)
summary(culeastm)
lowculeastm <- culeastm$ci[1,1]
upculeastm <- culeastm$ci[1,2]
coefculeastm <- culeastm$est[1]

culwestm <- RDestimate(cul ~ date | cntry, data = westm10, cutpoint = 0, bw = 14)
summary(culwestm)
lowculwestm <- culwestm$ci[1,1]
upculwestm <-culwestm$ci[1,2]
coefculwestm <- culwestm$est[1]



migeastr <- RDestimate(mig ~ date | cntry, data = eastr10, cutpoint = 0, bw = 14)
summary(migeastr)
lowmigeastr <- migeastr$ci[1,1]
upmigeastr <- migeastr$ci[1,2]
coefmigeastr <- migeastr$est[1]

migwestr <- RDestimate(mig ~ date | cntry, data = westr10, cutpoint = 0, bw = 14)
summary(migwestr)
lowmigwestr <- migwestr$ci[1,1]
upmigwestr <- migwestr$ci[1,2]
coefmigwestr <- migwestr$est[1]

ecoeastr <- RDestimate(eco ~ date | cntry, data = eastr10, cutpoint = 0, bw = 14)
summary(ecoeastr)
lowecoeastr <- ecoeastr$ci[1,1]
upecoeastr <- ecoeastr$ci[1,2]
coefecoeastr <- ecoeastr$est[1]

ecowestr <- RDestimate(eco ~ date | cntry, data = westr10, cutpoint = 0, bw = 14)
summary(ecowestr)
lowecowestr <- ecowestr$ci[1,1]
upecowestr <- ecowestr$ci[1,2]
coefecowestr <- ecowestr$est[1]

culeastr <- RDestimate(cul ~ date | cntry, data = eastr10, cutpoint = 0, bw = 14)
summary(culeastr)
lowculeastr <- culeastr$ci[1,1]
upculeastr <- culeastr$ci[1,2]
coefculeastr <- culeastr$est[1]

culwestr <- RDestimate(cul ~ date | cntry, data = westr10, cutpoint = 0, bw = 14)
summary(culwestr)
lowculwestr <- culwestr$ci[1,1]
upculwestr <-culwestr$ci[1,2]
coefculwestr <- culwestr$est[1]



hetcoefr2 <- c(coefmigeastl, coefmigwestl, coefecoeastl, 
              coefecowestl, coefculeastl, coefculwestl,
              coefmigeastm, coefmigwestm, coefecoeastm, 
              coefecowestm, coefculeastm, coefculwestm,
              coefmigeastr, coefmigwestr, coefecoeastr, 
              coefecowestr, coefculeastr, coefculwestr)
hetcoefr2 <- as.data.frame(hetcoefr2)

hetupr2 <- c(upmigeastl, upmigwestl, upecoeastl, upecowestl, 
            upculeastl, upculwestl,
            upmigeastm, upmigwestm, upecoeastm, upecowestm, 
            upculeastm, upculwestm,
            upmigeastr, upmigwestr, upecoeastr, upecowestr, 
            upculeastr, upculwestr)
hetupr2 <- as.data.frame(hetupr2)

hetlowr2 <- c(lowmigeastl, lowmigwestl, lowecoeastl, lowecowestl, 
             lowculeastl, lowculwestl,
             lowmigeastm, lowmigwestm, lowecoeastm, lowecowestm, 
             lowculeastm, lowculwestm,
             lowmigeastr, lowmigwestr, lowecoeastr, lowecowestr, 
             lowculeastr, lowculwestr)
hetlowr2 <- as.data.frame(hetlowr2)

hetplotr2 <- bind_cols(hetcoefr2,hetlowr2,hetupr2)
hetplotr2$varr <- as.character(c("Overall", "Overall", "Economy", "Economy", "Culture", "Culture",
                                 "Overall", "Overall", "Economy", "Economy", "Culture", "Culture",
                                 "Overall", "Overall", "Economy", "Economy", "Culture", "Culture"))
hetplotr2$both <- as.character(c("Eastern Left", "Western Left", "Eastern Left", 
                                 "Western Left", "Eastern Left", "Western Left",
                                 "Eastern Moderate", "Western Moderate", "Eastern Moderate", 
                                 "Western Moderate", "Eastern Moderate", "Western Moderate",
                                 "Eastern Right", "Western Right", "Eastern Right", 
                                 "Western Right", "Eastern Right", "Western Right"))
level_order <- c("Overall", "Economy", "Culture")

bothplot <- ggplot(hetplotr2, aes(x = varr, y = hetcoefr2, color = both)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 6, position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin = hetlowr2, ymax = hetupr2), 
                size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-3,3) +
  xlab("") +
  ylab("") +
  theme_light(base_size = 30) +
  scale_x_discrete(limits = level_order)
bothplot



# same without moderates, for better visibility of the relevant coefficients
# Define the coefficients, upper and lower limits for Left and Right values only
hetcoefr2 <- c(coefmigeastl, coefmigwestl, coefecoeastl, coefecowestl, coefculeastl, coefculwestl,
               coefmigeastr, coefmigwestr, coefecoeastr, coefecowestr, coefculeastr, coefculwestr)
hetcoefr2 <- as.data.frame(hetcoefr2)

hetupr2 <- c(upmigeastl, upmigwestl, upecoeastl, upecowestl, upculeastl, upculwestl,
             upmigeastr, upmigwestr, upecoeastr, upecowestr, upculeastr, upculwestr)
hetupr2 <- as.data.frame(hetupr2)

hetlowr2 <- c(lowmigeastl, lowmigwestl, lowecoeastl, lowecowestl, lowculeastl, lowculwestl,
              lowmigeastr, lowmigwestr, lowecoeastr, lowecowestr, lowculeastr, lowculwestr)
hetlowr2 <- as.data.frame(hetlowr2)

# Bind the data into a single dataframe
hetplotr2 <- bind_cols(hetcoefr2, hetlowr2, hetupr2)

# Add variable and group labels
hetplotr2$varr <- as.character(c("Overall", "Overall", "Economy", "Economy", "Culture", "Culture",
                                 "Overall", "Overall", "Economy", "Economy", "Culture", "Culture"))
hetplotr2$both <- as.character(c("Eastern Left", "Western Left", "Eastern Left", 
                                 "Western Left", "Eastern Left", "Western Left",
                                 "Eastern Right", "Western Right", "Eastern Right", 
                                 "Western Right", "Eastern Right", "Western Right"))

# Define the order for the x-axis
level_order <- c("Overall", "Economy", "Culture")


# Define custom colors for the groups
custom_colors <- c("Eastern Left" = "#CC0000", "Western Left" = "#FF6666", 
                   "Eastern Right" = "#004C99", "Western Right" = "#66B2FF")

# Create the plot
bothplot <- ggplot(hetplotr2, aes(x = varr, y = hetcoefr2, color = both)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 6, position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = hetlowr2, ymax = hetupr2), 
                size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-3, 3) +
  xlab("") +
  ylab("") +
  theme_light(base_size = 30) +
  scale_x_discrete(limits = level_order) +
  scale_color_manual(values = custom_colors)
bothplot





# respondents with migrant backgrounds
migsub <- ESS10 %>% 
  filter(brncntr == 2 | facntr == 2 | mocntr == 2)



rd_mig0 <- RDestimate(mig ~ date | cntry, data = migsub, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_mig0)
rd_eco0 <- RDestimate(eco ~ date | cntry, data = migsub, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_eco0)
rd_cul0 <- RDestimate(cul ~ date | cntry, data = migsub, cutpoint = 0, bw = 14, se.type = "HC1")
summary(rd_cul0)


